3.8. Radial PML#

by M. Wess, 2025#

This Notebook is part of the dualcellspaces documentation for the addon package implementing the Dual Cell method in NGSolve, as well as part of the td_evp documentation on the implementation of time-domain methods for resonance problems.

Derivation of the PML formulation#

We start from the time-domain Maxwell system in weak first-order \(E-H\) formulation

\[\begin{split} \partial_t \int_\Omega \tilde E\cdot \tilde e =\int_\Omega\mathrm{curl} \tilde H\cdot \tilde e,\\ \partial_t \int_\Omega \tilde H\cdot \tilde h = -\int_\Omega \mathrm{curl} \tilde E\cdot \tilde h \end{split}\]

for all test functions \(e,h\), where \(\Omega\) is some suitable domain.

Complex scaling#

To apply a PML we change the domain of integration to some complexified mapped domain \(\gamma(\Omega)\). Additionally we introduce the covariantly mapped fields

\[\begin{split} \begin{align} E &= J^{-T}\tilde E,& e &= J^{-T}\tilde e,\\ H &= J^{-T}\tilde H,& h &= J^{-T}\tilde h, \end{align} \end{split}\]

where \(J\) is the Jacobian matrix \(J=D\gamma\). After applying the transformation rule we obtain using

\[\begin{split} \mathrm{curl}\tilde E = \mathrm{curl}\left(J^{-T} E\right) = \frac{1}{\det J}J^T\mathrm{curl} E,\\ \mathrm{curl}\tilde H = \mathrm{curl}\left(J^{-T} H\right) = \frac{1}{\det J}J^T\mathrm{curl} H, \end{split}\]

the transformed weak formulation

\[\begin{split} \partial_t \int_\Omega J^{-T} E \cdot J^{-T} e \det J=\int_\Omega\mathrm{curl} H e,\\ \partial_t \int_\Omega J^{-T} H \cdot J^{-T} h \det J = -\int_\Omega \mathrm{curl} E h. \end{split}\]

Now the above formulation with a suitable complex coordinate transformation on a suitably truncated domain is our desired yields our desired PML formulation. However, as usual the complex transformation is derived in frequency domain with a suitable dependency on the frequency. Here, we skip the detour to frequency domain and indicate this dependency nonchalantly by a dependency on time-derivation \(\partial_t\).

Namely we choose a radial scaling

\[ \gamma(x) = x + d(\|x\|)\frac{x}{\|x\|}\frac{1}{\partial t}, \]

where \(d\) vanishes in the interior domain \(\Omega_{int}=B_R(0)\) for some \(R>0\) and is a suitably chosen damping profile in \(\Omega_{ext}=\Omega\setminus\overline{\Omega_{int}}=\overline{B_R(0)}^c\). A simple case is e.g., given by

\[\begin{split} d(r) = \begin{cases} 0,&r<R\\ \alpha (r-R),&r\geq R \end{cases} \end{split}\]

for a given damping constant \(\alpha>0\).

Thus, we have for the scaling Jacobian \(J\)

\[ \begin{align} J(x) &= I + \frac{d(\|x\|)}{\|x\|}I\frac{1}{\partial_t}+\left(\frac{d'(\|x\|)}{\|x\|^2}-\frac{d(\|x\|)}{\|x\|^3}\right)xx^T\frac{1}{\partial_t} \end{align} \]

We may rewrite the above using the (orthogonal) projections onto \(x\) and \(x^\perp\) using

\[ \begin{align} \Pi_{x} &= \frac{ xx^T}{\|x\|^2},&\Pi_x^\perp = I-\frac{xx^T}{\|x\|^2} \end{align} \]

to obtain

\[ J(x) = \left(1+d'(\|x\|)\frac{1}{\partial_t}\right)\Pi_x+\left(1+\frac{d(\|x\|)}{\|x\|}\frac{1} {\partial_t}\right)\Pi_x^\perp \]

and thus

\[ J^{-T}=\frac{\partial_t}{\partial_t+d'}\Pi_x+\frac{\partial_t \|x\|}{\|x\|\partial_t+d}\Pi_x^\perp \]
\[\begin{split} \begin{align} \partial_t J^{-1}J^{-T}\det J &= \frac{1}{\partial_t+d'}\frac{(\|x\|\partial_t+d)^2}{\|x\|^2}\Pi_x+(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t+\frac{2d\partial_t}{\|x\|}\frac{1}{\partial_t+d'}+\frac{d^2}{\|x\|^2}\frac{1}{\partial_t+d'}-\frac{d'\partial_t}{\partial_t+d'}\right)\Pi_x+(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t +\frac{2d}{\|x\|}-\frac{2dd'}{\|x\|}\frac{1}{\partial_t+d'} +\frac{d^2}{\|x\|^2}\frac{1}{\partial_t+d'} -d'+\frac{d'd'}{\partial_t+d'} \right)\Pi_x +(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t-d'+\frac{2d}{\|x\|} +\left(-\frac{2dd'}{\|x\|} +\frac{d^2}{\|x\|^2} +d'd' \right)\frac{1}{\partial_t+d'} \right)\Pi_x +(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t-d'+\frac{2d}{\|x\|} +\frac{\left(d'-\frac{d}{\|x\|} \right)^2}{\partial_t+d'} \right)\Pi_x +(\partial_t+d')\Pi_x^\perp \end{align} \end{split}\]

The full system#

introducing the auxiliary unknowns

\[E_1:= \frac{d'-\frac{d}{\|x\|}}{\partial_t+d'} E\]

and similarly for \(H_1\) leads to

\[\begin{split} \partial_t E_1 = (d'-\frac{d}{\|x\|})E-d'E_1\\ \partial_t H_1 = (d'-\frac{d}{\|x\|})H-d'H_1. \end{split}\]

Thus, we obtain the full PML system

\[\begin{split} \begin{align} \partial_t \int_\Omega E\cdot e &=\int_\Omega(d'-\frac{2d}{\|x\|})\Pi_x E\cdot e+\int_\Omega(-d'+\frac{d}{\|x\|})\Pi_x E_1\cdot e-\int_\Omega d' \Pi_x^\perp E\cdot e+\int_\Omega\mathrm{curl} H e,\\ \partial_t\int_\Omega E_1 \cdot e_1 &= \int (d'-\frac{d}{\|x\|})E\cdot e_1-\int d' E_1\cdot e_1\\ \partial_t \int_\Omega H\cdot h &= \int_\Omega(d'-\frac{2d}{\|x\|})\Pi_x H\cdot h+\int_\Omega(-d'+\frac{d}{\|x\|})\Pi_x H_1\cdot h-\int_\Omega d' \Pi_x^\perp H\cdot h-\int_\Omega \mathrm{curl} E h\\ \partial_t\int_\Omega H_1\cdot h_1 &= \int (d'-\frac{d}{\|x\|})H\cdot h_1-\int d' H_1\cdot h_1 \end{align} \end{split}\]

Implementation#

We start by doing the necessary imports

from ngsolve import *
import dualcellspaces as dcs
from ngsolve.webgui import Draw
from netgen.occ import *
import numpy as np
import matplotlib.pyplot as pl
R = 0.8     # inner sphere radius
R_pml = 1.1     # outer sphere radius

maxh = 0.06

order = 0

alpha = 10


# Create outer and inner spheres (both centered at origin)
outer_sphere = Sphere(Pnt(0, 0, 0), R_pml)
inner = Sphere(Pnt(0, 0, 0), R)

pml = outer_sphere-inner
pml.name = "pml"
inner.name = "inner"

geo = OCCGeometry(Glue([pml,inner]))
mesh = Mesh(geo.GenerateMesh(maxh = maxh))

#mesh.Curve(order)
Draw(mesh)
print(mesh.GetMaterials())
('pml', 'inner')
fes_E = dcs.HCurlDualCells3D(mesh, order=order)
fes_H = dcs.HCurlPrimalCells3D(mesh,order=order)

print("total DoFs:",fes_E.ndof+fes_H.ndof)
gfE, gfH = GridFunction(fes_E), GridFunction(fes_H)
gfE1, gfH1 = GridFunction(fes_E), GridFunction(fes_H)


#integral symbols with special integration rules
dxE = dx(intrules=fes_E.GetIntegrationRules())
dxH = dx(intrules=fes_H.GetIntegrationRules())

dxw = dx(intrules=dcs.GetIntegrationRules(2*order+4))
dSw = dx(element_boundary=True,intrules=dcs.GetIntegrationRules(2*order+4))


#mixed bilinear form
E,dE = fes_E.TnT()
H,dH = fes_H.TnT()

normal = specialcf.normal(3)
bf_mixed = BilinearForm(E*curl(dH)*dxw+E*Cross(dH,normal)*dSw, geom_free=True).Assemble().mat
total DoFs: 489292
pmldofs_E = fes_E.GetDofs(mesh.Materials('pml'))
pmldofs_H = fes_H.GetDofs(mesh.Materials('pml'))
# coefficients for mass matrices
vx = CF((x,y,z))
r = Norm(vx)

Pi = OuterProduct(vx/r,vx/r)
Pi_perp = Id(3)-Pi


d = mesh.MaterialCF({"inner" : 0, "pml" : R+(r-R)*alpha})

d_ = d.Diff(r)


#Draw(d,mesh,draw_surf=False, clipping={"y":0, "z":-1},
#      settings = {"Objects" : {"Clipping Plane" : True}}, euler_angles=[0,40,0])

#Draw(d_,mesh,draw_surf=False, clipping={"y":0, "z":-1},
#      settings = {"Objects" : {"Clipping Plane" : True}}, euler_angles=[0,40,0])
#prepare mass operators
with TaskManager():
    # (inverse) mass
    massH_inv = fes_H.Mass(1).Inverse()
    massE_inv = fes_E.Mass(1).Inverse()
    
    # inverse mass in pml only   
    massH_inv_pml = fes_H.Mass(1).Inverse(freedofs = pmldofs_H)
    massE_inv_pml = fes_E.Mass(1).Inverse(freedofs = pmldofs_E)    
    
    # mass matrices for pml damping
    dampEe = BilinearForm((d_-2*d/r)*(Pi*E)*dE*dxE-d_*Pi_perp*E*dE*dxE).Assemble().mat
    dampE1e = BilinearForm((-d_+d/r)*(Pi*E)*dE*dxE).Assemble().mat
    dampEe1 = BilinearForm((d_-d/r)*E*dE*dxE).Assemble().mat
    dampE1e1 = BilinearForm(-d_*E*dE*dxE).Assemble().mat

    dampHh = BilinearForm((d_-2*d/r)*(Pi*H)*dH*dxH-d_*Pi_perp*H*dH*dxH).Assemble().mat
    dampH1h = BilinearForm((-d_+d/r)*(Pi*H)*dH*dxH).Assemble().mat
    dampHh1 = BilinearForm((d_-d/r)*H*dH*dxH).Assemble().mat
    dampH1h1 = BilinearForm(-d_*H*dH*dxH).Assemble().mat
    
def estimate_tau(mat, maxsteps = 2000, tol = 1e-9):   
    vec = mat.CreateColVector()
    vec.SetRandom()
    vec*=1/vec.Norm()
    tmp = vec.CreateVector()
    lam = 0
    
    for i in range(maxsteps):
        #print(i,end='\r')
        tmp.data = mat * vec
        
        lamnew = InnerProduct(tmp,vec)
        tau = 2/sqrt(lamnew)
        n = tmp.Norm()
        res = (tmp-lam*vec).Norm()/lamnew/n
        #print(res)
        tmp *= 1/n
        if res<tol: return tau
        vec.data = tmp
        lam = lamnew
    print("did not converge, last res = ",res)
    return tau

tau = estimate_tau(massE_inv@bf_mixed.T@massH_inv@bf_mixed)
tau *= 0.9
print("estimated tau = ", tau)
estimated tau =  0.010683411757906781
#set initial data

a = exp(-r**2*20)


E0 = CF((a.Diff(y),-a.Diff(x),0))  #E0 = curl(a*(0,0,1))
rhs = LinearForm(E0*dE*dxE).Assemble().vec
gfE.vec.data = massE_inv*rhs

#Draw solution
#scenee = Draw (gfE.Operator("altshape"), mesh, "E", order=2, draw_surf=False, clipping={"y":0, "z":-1},
#               min=-0.5,max=0.5,autoscale=False,points=dcs.GetWebGuiPoints(2), vectors=True,
#               settings = {"eval": 0, "Objects" : {"Clipping Plane" : True}}, euler_angles=[0,30,0])

# or save solution in multivector for animation
#gfE_anim = GridFunction(fes_E,multidim=0)

from time import time

drawevery = 70
tend = 1.1

t = 0.
i = 0

gfH.vec[:]=0
#gfE.vec[:]=0

gfH1.vec[:]=0
gfE1.vec[:]=0

tmpE = gfE.vec.CreateVector()
tmpH = gfH.vec.CreateVector()

tmpE_ = gfE.vec.CreateVector()
tmpH_ = gfH.vec.CreateVector()



now = time()
timepassed = 0


damp = True

with TaskManager():
    gfH.vec.data += -tau/2*massH_inv@bf_mixed*gfE.vec     #initial data in pml is zero, so first step does not need damping terms
    while t<tend:
        if i%drawevery == 0:
            timepassed += time()-now
            #scenee.Redraw() # for immediate drawing
            #gfE_anim.AddMultiDimComponent(gfE.vec) #for saving to multidim function
            Draw (gfE.Operator("altshape"), mesh, "E0", order=1, draw_surf=False, clipping={"y":0, "z":-1},
               min=-1,max=1,autoscale=False,points=dcs.GetWebGuiPoints(1), 
               settings = {"eval": 0, "Objects" : {"Clipping Plane" : True}}, euler_angles=[0,30,0])            
            print(" time = {}, norm = {}, step = {}, {:e} dofs/s".format(t,gfE.vec.Norm(),i,i*(fes_E.ndof+fes_H.ndof)/timepassed),end="")

            now = time()
        t+=tau       
        i+=1

        if damp:
            tmpE.data = dampE1e1*gfE1.vec+dampEe1*gfE.vec
            tmpE_.data = bf_mixed.T*gfH.vec+dampEe*gfE.vec+dampE1e*gfE1.vec
            gfE.vec.data += tau*massE_inv*tmpE_
            gfE1.vec.data += tau*massE_inv_pml*tmpE


            tmpH.data = dampH1h1*gfH1.vec+dampHh1*gfH.vec
            tmpH_.data = -bf_mixed*gfE.vec+dampHh*gfH.vec+dampH1h*gfH1.vec
            gfH.vec.data += tau*massH_inv*tmpH_
            gfH1.vec.data += tau*massH_inv_pml*tmpH
        
        else:
            gfE.vec.data += tau*massE_inv@bf_mixed.T*gfH.vec
            gfH.vec.data += -tau*massH_inv@bf_mixed*gfE.vec


        
 time = 0.0, norm = 4.573211163926184, step = 0, 0.000000e+00 dofs/s
 time = 0.747838823053474, norm = 2.750803648279417, step = 70, 1.382374e+07 dofs/s
#print("finshed with an average of {:e} dofs/s".format(i*(fes_E.ndof+fes_H.ndof)/timepassed),end="")

#Draw (gfE_anim, mesh, "E0", order=2, draw_surf=False, clipping={"y":0, "z":-1},
#               min=-1,max=1,autoscale=False,points=dcs.GetWebGuiPoints(1), 
#               animate = True,
#               settings = {"eval": 0, "Objects" : {"Clipping Plane" : True}}, euler_angles=[0,30,0])